Fit a Generalized Linear Model for Peak Offset
clearvars -except T Gr
if exist('Gr','var')==0 % If table not in workspace, load it.
Gr = getfield(load(defaults.files('single_channel_stats_mlx_table')),'Gr');
end
We would like to fit a model that predicts the Peak Offset parameter that we extracted in analyze.stat.fit_gauspuls. We are primarily interested in the fixed effect of:
- GroupID : Is there an effect of whether the rat has an ischemia or not?
- PostOpDay : Is there a trend that correlates with recovery over time? Rats performed the task daily between post-op days 3 to 28.
- Area : Either "RFA" or "CFA"
We also think that the following random effects could have an impact on our model intercept:
- ChannelID : We expect some channels intrinsically to have an earlier or later peak relative to others.
- AnimalID : We expect some animals to have a different aggregate distribution of channels.
Finally, at the Channel level, we also think there could be an influence of coefficient for the following terms:
- Error_SS : The residual sum-of-squares from individual gauspuls trial fits.
- Duration : The total time that it took for the reach-to-grasp behavior to finish execution.
Fit model for PeakOffset
We can fit a generalized linear mixed-effects model (GLME) using the Matlab fitglme function that requires the 'Statistics and Machine Learning Toolbox', version 12.0. This function takes a model specified in Wilkinson Notation, which defines Fixed and Random effects. We will fit using the ApproximateLaplace method, which allows us to compare model likelihoods if needed (as opposed to a pseudolikelihood-based method).
glme_Offset = fitglme(Gr,"PeakOffset~GroupID*PostOpDay*Area+(1+Error_SS+Duration|ChannelID)+(1|AnimalID)", ...
'FitMethod','ApproximateLaplace');
disp(glme_Offset);
Model Fixed-Effect Coefficients
All probabilities are computed using two-tailed t-tests for a distribution with an estimated degrees of freedom of 20,213.
- (Intercept) : The term with the greatest effect on the model is the average offset value. On average, peak times precede the Grasp by approximately 150-ms (coefficient estimate of -0.15056 seconds). This effect is statistically significant, with tStat=-11.327 and p=1.19x10^-29. For other coefficients in this list, the tStat value is omitted, but you can check the table above in the corresponding column.
- GroupID[Intact] : All "dummy" variables that are created for fixed-effects involving a categorical variable have n-1 coefficients estimated (since they use the "reference" method). In this case, Ischemia is set to zero (the reference) and the coefficient for Intact is estimated. The coefficient estimate for the effect is very small, and the p-value (p=0.77) indicates that this is not a significant model term.
- Area[CFA] : As with GroupID, Area only has the coefficient for one level computed- CFA (the intact hemisphere primary motor cortex homolog, which is ipsilateral to the reaching forelimb). However, unlike GroupID, there is a significant effect for this main effect coefficient (p=0.00075), which indicates that channels in CFA have a peak offset that on average is approximately 60-ms (0.0613 seconds) closer to the Grasp than those in RFA. This makes sense: we expect RFA, a putative premotor homolog, to in general have peak-excitation earlier in time relative to the grasp behavior.
- PostOpDay : This is a continuous fixed-effect, and the model indicates it is also a significant term (p=0.0051). The coefficient value of 0.0022 (i.e. 2-ms per day) seems a little bit low, but over the course of 30 days that can be a 60-ms shift, which is substantial.
- GroupID[Intact]*Area[CFA], GroupID[Intact]*PostOpDay : Neither term tests statistically significant (p=0.38, and p=0.83 respectively). These terms indicate that the change in timing between CFA-RFA channels, as well as the slope of expected peak offset by day, does not change based on whether the rat has an ischemic infarct or not.
- Area[CFA]*PostOpDay : This interaction term tests significant at alpha=0.05 level, with a negative coefficient indicating that for CFA channels, there is a decline of approximately 2-ms per day that seems to match the expected overall fixed-effect of PostOpDay. This indicates that we expect the peak offset time to stay relatively constant in CFA, whereas in RFA (for either group), we expect the value to shift closer to the Grasp.
- Area[CFA]*GroupID[Intact]*PostOpDay : The three-way interaction is not significant, meaning there is no deviation in slope of expected Peak Offset by post-op day that is specific to one combination of Group and Area.
Check model fit based on explained variance
disp(glme_Offset.Rsquared);
One concern after checking the Rsquared values is that our model doesn't really capture a lot of the data.
Check model fit using scatter plots
analyze.stat.scatter_var(Gr,'PeakOffset_ms','AddLabel',false,'YLab',"Peak-Offset");
We can see from the scatter that the reason our model fit is poor is because the range of peak times relative to behavior on any given trial is quite large relative to any overall trends in average peak times either by group, area, or post-op day. In general, there seems to be a robust difference between RFA and CFA peaks, with RFA peaks tending to be earlier in general than CFA peaks, which is what we expected.
Plot residuals to see if we should restrict range
analyze.stat.panelized_residuals(glme_Offset);
Based on the distribution of residuals, we see that there is a multi-phasic histogram of residuals (which isn't good), and that a lot of our fitted values are on the range [-0.2 to 0.1] approximately.
Set the included Peak Offset cutoff points
We can look at the distribution of PeakOffset so we have an idea of what our response variable looks like in general (should have done this first).
figure('Name','Peak Offset Observed Distribution','Color','w','Units','Normalized','Position',[0.2 0.2 0.35 0.35]);
histogram(Gr.PeakOffset,100,'DisplayName','Counts');
xlabel('Peak Offset (sec)','FontName','Arial','Color','k','FontSize',14,'FontWeight','bold');
ylabel('Count','FontName','Arial','Color','k','FontSize',14,'FontWeight','bold');
title('Total Distribution of Peak Offset Times','FontName','Arial','FontSize',16,'FontWeight','bold','Color','k');
Add lines to the histogram to denote where we initialize our "cutoff" for exclusions in next section. These should move as sliders are adjusted.
yy = [0 1000];
x1 = [-1.00 -1.00];
x2 = [ 0.75 0.75];
h_min = line(x1,yy,'Color',[1.0 0.0 0.0],'LineWidth',2,'LineStyle',':','DisplayName','min\_peak\_offset');
h_max = line(x2,yy,'Color',[0.8 0.2 0.2],'LineWidth',2,'LineStyle',':','DisplayName','max\_peak\_offset');
legend(gca,'Location','Best');
Set min_peak_offset and max_peak_offset by dragging the sliders to the desired range.
min_peak_offset = -0.8;
max_peak_offset = 0.4;
set(h_min,'XData',ones(1,2).*min_peak_offset);
set(h_max,'XData',ones(1,2).*max_peak_offset);
drawnow;
Fit the model with exclusions applied
exclusions = (Gr.PeakOffset<min_peak_offset) | (Gr.PeakOffset>max_peak_offset);
glme_Offset_subset_times = fitglme(Gr,"PeakOffset~GroupID*PostOpDay*Area+(1+Error_SS+Duration|ChannelID)+(1|AnimalID)", ...
'FitMethod','ApproximateLaplace','Exclude',exclusions);
disp(glme_Offset_subset_times);
After applying the exclusions, our model has roughly the same parameter explanations.
disp(glme_Offset_subset_times.Rsquared);
And actually the fit is worse than before.
Look at the included data to make sure we are trying to fit a reasonable manifold
Depending upon how the data exclusion/inclusion criteria were applied, we could end up with an ill-posed regression.
analyze.stat.scatter_var(Gr(~exclusions,:),'PeakOffset_ms',...
'AddLabel',false,...
'YLab',"Peak-Offset",...
'ShowChannels',true, ... % Plot individual channels
'ScatterType','mean',... % Look at channel means
'ShowCB',true,... % Do look at individual channel confidence bands
'YLim',[min_peak_offset max_peak_offset].*1e3);
From this plot, we should be able to tell pretty quickly if we have been too restrictive in setting the mean boundaries: the total number of samples will seem sparse in one of the panels. The optional parameter specifications in the call to analyze.stat.scatter_var alter what we are plotting a little bit: instead of individual trials, data points represent Channel Averages across all trials from a given recording session. Confidence bands show 95% bounds on Animal regression estimates (regression using all channels from a single animal to estimate expected (average) Peak-Offset on a given Post-Op Day). Unfortunately our estimates for Intact animals are all over the place, which might be why the model fits in general are pretty noisy or hard to interpret.
Re-plot Residuals for Peak-Offset fit using exclusions
analyze.stat.panelized_residuals(glme_Offset_subset_times);
It's still not very good. We can try to restrict to only sites that are Forelimb-Related or non-Forelimb-Related to see if it helps.
Restrict to only Forelimb-related
Using intracortical microstimulation (ICMS), we approximately co-registered the somatic territory neareste to a given microwire shank. Each channel therefore has a categorical tag indicating the corresponding movement representation: Distal Forelimb (DF); Proximal Forelimb (PF); Forelimb Boundary (DF-PF); Other (Neck, Vibrissae, Trunk, Shoulder, Mouth; "O"); or No Response (NR).
We can create two categories: one for forelimb-related and one for "other" and then make models grouped by those categories.
Create ICMS_Type variable and only use "Forelimb" sites
% Create new variable and use it for exclusions
forelimb_channels = (Gr.ICMS=="DF") | (Gr.ICMS=="PF") | (Gr.ICMS=="DF-PF");
Gr.ICMS_Type = categorical(forelimb_channels,[false true],["Other", "Forelimb"]);
exclusions = Gr.ICMS_Type=="Other";
Run GLME fit, including full range of Peak Offset times
glme_Offset_forelimb = fitglme(Gr,"PeakOffset~GroupID*PostOpDay*Area+(1+Error_SS+Duration|ChannelID)+(1|AnimalID)", ...
'FitMethod','ApproximateLaplace','Exclude',exclusions);
disp(glme_Offset_forelimb);
disp(glme_Offset_forelimb.Rsquared);
The model explains these data a little worse than the very first model. Interestingly, when restricted to only forelimb sites, no model fixed effect is significant at the alpha=0.05 level except for the (Intercept) term.
Visualize scatter of the data used in model
analyze.stat.scatter_var(Gr(~exclusions,:),'PeakOffset_ms','AddLabel',false,'YLab',"Peak-Offset");
Again, grey bands indicate 95% confidence bands for all trial data and each scatter point is a single trial value for a single channel. Nothing too interesting here.
Plot residuals of Forelimb-only Peak-Offset model
analyze.stat.panelized_residuals(glme_Offset_forelimb);
And again, the residuals are less than ideal.
Fit second model based on only non-forelimb sites
It may be the case that plasticity occurs via recruitment of non-forelimb sites, in which case those may actually be the more interesting channels to look at.
exclusions = Gr.ICMS_Type=="Forelimb";
glme_Offset_other = fitglme(Gr,"PeakOffset~GroupID*PostOpDay*Area+(1+Error_SS+Duration|ChannelID)+(1|AnimalID)", ...
'FitMethod','ApproximateLaplace','Exclude',exclusions);
disp(glme_Offset_other);
The model terms generally agree with the original model.
disp(glme_Offset_other.Rsquared);
This model actually fits the data better than the original model.
Verify that the included data is appropriate
analyze.stat.scatter_var(Gr(~exclusions,:),'PeakOffset_ms','AddLabel',false,'YLab',"Peak-Offset");
The scatter plot shows us that we just don't have enough data to look at the model this way, since there are very few non-forelimb sites in the Intact CFA grouping. So it would be problematic to make too much of a conclusion from this restricted grouping.
analyze.stat.panelized_residuals(glme_Offset_other);
Fit Generalized Linear Model for Envelope Bandwidth
We use the same structure as for fitting Peak Offset, but this time we will use a log link since we are dealing with frequencies.
glme_BW = fitglme(Gr,"EnvelopeBW~GroupID*Area*PostOpDay+(1+Error_SS+Duration|ChannelID)+(1|AnimalID)",...
'link','log','FitMethod','ApproximateLaplace');
disp(glme_BW);
disp(glme_BW.Rsquared);
This model suggests there is a significant effect of Group on Envelope Bandwidth; however, when we go to look at the effect size it really is quite small relative to what we could think of making a conceivable difference in terms of a frequency parameter. It's possible that a change on the order of tenths of a Hz relates to some physiological process of interest, so we will plot these data and see if we can get a better intuition of what the parameter significance means.
analyze.stat.scatter_var(Gr,'EnvelopeBW',...
'AddLabel',false,...
'YLab',"Envelope Bandwidth",...
'formatspec','%6.3f',...
'ShowChannels',true, ... % Plot individual channels
'ScatterType','mean',... % Look at channel means
'ShowCB',true,... % Do look at individual channel confidence bands
'YLim',[0 2]);
We can definitely see that there is some kind of by-day tendency across the board for a (small) increase in Envelope Bandwidth. A change from 0.8-Hz to 0.9-Hz does seem pretty small, but the difference in period between frequencies at that scale is actually substantial, so even that small change would correspond to a difference of approximately 100-ms. Therefore, it may make sense to do a transformation on this response variable and look at it in terms of the period to see if the model is better fit.
Create transformed response variable and fit response
First, create a variable called Envelope_Period for the purpose of visualizing on scatter plots, in units of milliseconds. Then, create a normalized log-transformed version of that variable, which should be better for fitting a statistical model.
Gr.Envelope_Period = 1000./Gr.EnvelopeBW;
Gr.Properties.VariableUnits{28} = 'ms';
Gr.Envelope_Period_Transformed = (log(Gr.Envelope_Period) - nanmean(log(Gr.Envelope_Period)))./nanstd(log(Gr.Envelope_Period));
Gr.Properties.VariableUnits{29} = 'z-score';
figure('Name','Difference in Distribution from Transformation');
subplot(1,2,1); histogram(Gr.Envelope_Period); title('Period (ms)'); subplot(1,2,2); histogram(Gr.Envelope_Period_Transformed); title('Transformed Distribution');
Whereas the original distribution is right-skewed, the new distribution is approximately "balanced" (although a kstest indicates that it is not statistically "normal").
Fit generalized linear model of transformed period variable
Note: since we are working with the transformed variable, we will not use the log link function.
glme_T_env = fitglme(Gr,"Envelope_Period_Transformed~GroupID*Area*PostOpDay+(1+Error_SS+Duration|ChannelID)+(1|AnimalID)",...
'FitMethod','ApproximateLaplace');
disp(glme_T_env);
The main fixed-effects are significant, but none of the interaction effects test as significant.
disp(glme_T_env.Rsquared);
When we plot the scatter, we will use the non-transformed version so that the data values make more sense.
analyze.stat.scatter_var(Gr,...
'Envelope_Period',...
'AddLabel',false,...
'YLab',"Envelope Period",...
'ShowChannels',true, ... % Plot individual channels
'ScatterType','mean',... % Look at channel means
'ShowCB',true,...
'YLim',[0 2500]); % Even a value of 2500-ms is long on the timescale we are interested in
This trend pretty much matches what we describe above. It may relate to the Duration random effect: if the behavior in general gets faster over time, this could be due to some common-source. We can look at statistics for the Duration random effect:
[~,~,STATS] = randomEffects(glme_T_env);
K = [STATS.Name] == "Duration";
pVal = coefTest(glme_T_env,zeros(1,8),0,'REContrast',double(K'));
disp(pVal);
The coefficients for Duration trend towards significance, suggesting that this effect could relate to what we have observed.
Fit generalized linear model for Center Frequency
We use the same structure as for fitting Peak Offset, but this time we will use a log link since we are dealing with frequencies.
glme_FC = fitglme(Gr,"PeakFreq~GroupID*Area*PostOpDay+(1+Error_SS+Duration|ChannelID)+(1|AnimalID)",...
'link','log','FitMethod','ApproximateLaplace');
disp(glme_FC);
disp(glme_FC.Rsquared);
analyze.stat.scatter_var(Gr,'PeakFreq',...
'AddLabel',false,...
'formatspec','%6.3f',...
'YLab','Center Frequency',...
'ShowChannels',true, ... % Plot individual channels
'ScatterType','mean',... % Look at channel means
'ShowCB',true,... % Do look at individual channel confidence bands
'YLim',[0 2]);
analyze.stat.panelized_residuals(glme_FC);
